% SABRE experiment simulation for Eibe Duecker and Christian Griesinger.
% Set to reproduce Figure 3b from 
%
%        http://dx.doi.org/10.1021/ja903601p
%
% i.kuprov@soton.ac.uk

function sabre_pyridine()

% Spin system
sys.isotopes={'1H','1H','1H','1H','1H','1H','1H'};

% Chemical shifts
inter.zeeman.scalar={8.54 7.44 7.86 7.44 8.54 -23.5 -23.5};

% Couplings inside pyridine
inter.coupling.scalar=cell(7);
inter.coupling.scalar{1,2}= 4.88; 
inter.coupling.scalar{4,5}= 4.88;
inter.coupling.scalar{1,4}= 1.00; 
inter.coupling.scalar{2,5}= 1.00;
inter.coupling.scalar{1,3}= 1.84; 
inter.coupling.scalar{3,5}= 1.84;
inter.coupling.scalar{1,5}=-0.13;
inter.coupling.scalar{2,3}= 7.67; 
inter.coupling.scalar{3,4}= 7.67;
inter.coupling.scalar{2,4}= 1.37;

% Couplings of the hydride group - please check
inter.coupling.scalar{1,6}=1.12;
inter.coupling.scalar{1,7}=1.02;
inter.coupling.scalar{6,7}=7.00;

% Magnetic fields
polarization_field=25e-3;
nmr_field=7.05;

% Basis set
bas.mode='complete';

% Do the housekeeping
sys.magnet=polarization_field;
spin_system=create(sys,inter);
spin_system=basis(spin_system,bas);

% Get the Hamiltonian superoperator
H=h_superop(assume(spin_system,'nmr'));

% Start in a singlet state
rho=singlet(spin_system,6,7);

% Evolve the system for 2.5 seconds
[timestep,nsteps]=stepsize(H,2.5);
rho=evolution(spin_system,H,[],rho,timestep,nsteps,'final');

% Disconnect the parahydrogen
[H,rho]=decouple(spin_system,H,rho,[6 7]);

% Determine the optimal time step
[timestep,nsteps]=stepsize(H,2.5);

% Evolve the system for a further 2.5 seconds
rho=evolution(spin_system,H,[],rho,timestep,nsteps,'final');

% Split the Hamiltonian superoperator
Hz=h_superop(assume(spin_system,'nmr','zeeman'));
Hc=h_superop(assume(spin_system,'nmr','couplings'));
Hz=decouple(spin_system,Hz,[],[6 7])/polarization_field;
Hc=decouple(spin_system,Hc,[],[6 7]);

% Lift the field exponentially in 1024 steps over 5 seconds
 for field=exp(linspace(log(polarization_field),log(nmr_field),1024));
     rho=step(spin_system,Hc+field*Hz,rho,5/1024);
 end

% Set the field to NMR field
H=Hc+nmr_field*Hz;
 
% Evolve the system for a further second in high field
[timestep,nsteps]=stepsize(H,1.0);
rho=evolution(spin_system,H,[],rho,timestep,nsteps,'final');

% Set NMR experiment parameters
parameters.offset=0;
parameters.spins={'1H'};
parameters.sweep=2000;
parameters.npoints=8196;
parameters.zerofill=65536;
parameters.axis_units='Hz';
parameters.invert_axis=1;

% Set the detection state
coil=state(spin_system,'L+',parameters.spins{1});
coil=coil/norm(coil);

% Get the pulse operator
Lp=operator(spin_system,'L+',parameters.spins{1});

% Apply pi/2 pulse on Y axis
rho=step(spin_system,(Lp-Lp')/2i,rho,pi/2);

% Detect the magnetization
fid=evolution(spin_system,H,coil,rho,1/parameters.sweep,...
              parameters.npoints-1,'observable');

% Apodization
fid=apodization(fid,'exponential-1d',20);

% Fourier transform
spectrum=fftshift(fft(fid,parameters.zerofill));

% Plotting
plot_1d(spin_system,real(spectrum),parameters);

end

